Mapping underserved communities

Below, we have provided instructions for replicating our methodology for mapping underserved communities in Tampa Bay. To view instructions for downloading the necessary source data from EJScreen, see Getting demographic data.

Load the required R packages (install first as needed).

library(sf)
library(mapview)
library(dplyr)
library(RColorBrewer)

Load and view the data.

load(file = 'data/dattb.RData')
mapview(dattb)

You can see that some census tracts are only representative of the bay. We can clean up this data by retaining only the tracts in which the total population recorded in the latest American Community Survey (“ACSTOTPOP”) was above zero. This will remove tracts in which no people reside (e.g., large waterbodies, parks, or other natural areas).

In line with EPA recommendations, we will use the following five demographic variables to identify underserved communities:

• Percent classified as low income

• Percent unemployed

• Percent linguistically isolated (i.e., limited English speaking)

• Percent with less than a high school education

• Percent with low life expectancy

The EPA recommends flagging communities that fall within the 80th percentile (or higher) nationally as potentially disadvantaged. However, meeting this threshold in one variable alone is not necessarily an appropriate measure of disadvantage (e.g., a community of predominantly wealthy retirees would meet the unemployment threshold). For TBEP’s Equity Strategy, we define an underserved community as one that meets at least two of these thresholds, which we believe reduces the number of errors in identification while recognizing that a community does not need to meet the threshold of every indicator to face significant challenges.

Run the code below to (1) remove unpopulated census tracts, (2) count the number of demographic thresholds met in each tract, and (3) identify which tracts will be classified as “underserved” communities.

dattbindex <- dattb %>%
  filter(ACSTOTPOP > 0) %>%
  mutate(threshold_income = ifelse(P_LWINCPCT >= 80, 1, 0),
         threshold_unempl = ifelse(P_UNEMPPCT >= 80, 1, 0),
         threshold_lingui = ifelse(P_LNGISPCT >= 80, 1, 0),
         threshold_educat = ifelse(P_LESHSPCT >= 80, 1, 0),
         threshold_lifexp = ifelse(P_LIFEEXPCT >= 80, 1, 0)) %>%
  rowwise() %>%
  select(matches('^threshold|^ID')) %>% 
  mutate(threshold_N = sum(threshold_income,threshold_unempl,threshold_lingui,threshold_educat,threshold_lifexp, na.rm = TRUE)) %>%
  mutate(underserved = ifelse(threshold_N > 1, "Yes", "No"))

View the first five rows to see how the calculations have played out.

head(dattbindex)

Simple feature collection with 6 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -82.45127 ymin: 28.03256 xmax: -82.39752 ymax: 28.05464
Geodetic CRS:  WGS 84
# A tibble: 6 × 9
# Rowwise: 
  ID          thresh…¹ thres…² thres…³ thres…⁴ thres…⁵                     Shape
  <chr>          <dbl>   <dbl>   <dbl>   <dbl>   <dbl>             <POLYGON [°]>
1 12057000101        1       0       0       0      NA ((-82.42619 28.04506, -8…
2 12057000102        1       1       1       1       0 ((-82.42633 28.03258, -8…
3 12057000201        1       1       1       1       1 ((-82.4512 28.04362, -82…
4 12057000202        1       1       1       1      NA ((-82.44289 28.05294, -8…
5 12057000301        1       0       1       1       1 ((-82.45127 28.03364, -8…
6 12057000302        1       0       0       1       1 ((-82.43595 28.0391, -82…
# … with 2 more variables: threshold_N <dbl>, underserved <chr>, and
#   abbreviated variable names ¹​threshold_income, ²​threshold_unempl,
#   ³​threshold_lingui, ⁴​threshold_educat, ⁵​threshold_lifexp

View by Census Tract

View a map showing the number of thresholds met per census tract (you may adapt to a different color scale of your choice).

mapview(dattbindex, zcol = "threshold_N", col.regions = brewer.pal(6, "Reds"), layer.name = "No. of Thresholds Met")

View the tracts that meet our definition of underserved communities. The areas in red are those that rank in the 80th percentile (or greater) nationally in 2 or more of the demographic screening variables. They will serve as priority areas for increasing the equitable distribution of benefits from TBEP’s environmental programs.

mapview(dattbindex, zcol = "underserved", col.regions = list("gray","red"), layer.name = "Underserved Communities")

You can save this final data as an RData object for future use.

# save the layer as an RData object
save(dattbindex, file = 'data/dattbindex.RData')

View by Drainage Basin

While the census tract delineations provide a practical map for identifying target neighborhoods and stakeholders that may be working for or in these underserved communities, it is not well-aligned with geographic delineations that are most informative for planning conservation and and restoration projects across the watershed.

For environmental planning purposes, TBEP will thus characterize unique water body assessment units according to the proportion of each unit containing underserved census tracts. This approach serves as a bridge between social and ecological units relevant to different stakeholders across Tampa Bay.

The Florida Department of Environmental Protection provides the Waterbody IDs (WBIDs) dataset, available here, which includes polygons delineating the drainage basins surrounding water body assessment units.

Load the drainage basin data.

dbasins <- st_read('https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/WBIDS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson')

Reading layer `OGRGeoJSON' from data source 
  `https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/WBIDS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson' 
  using driver `GeoJSON'
Simple feature collection with 6788 features and 17 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.6349 ymin: 24.39631 xmax: -79.97431 ymax: 31.00065
Geodetic CRS:  WGS 84

Intersect the layer with the Tampa Bay watershed.

load(file = 'data/tbshed.RData')

sf_use_s2(FALSE)

dattbdbasins <- dbasins %>% 
  st_intersection(tbshed)

Since we will be doing area calculations, we need to set a projected coordinate system for the drainage basin and census tract layers. The most appropriate CRS for Tampa Bay is EPSG 6443.

dattbdbasins <- st_transform(dattbdbasins, crs = 6443)
dattbindex <- st_transform(dattbindex, crs = 6443)

Keep only the land area within each drainage basin and census tract. We can exclude water areas by intersecting the layers with Florida’s shoreline (available here from the Florida Fish and Wildlife Conservation Commission).

flshore <- st_read('https://atoll.floridamarine.org/arcgis/rest/services/FWC_GIS/OpenData_Shoreline/MapServer/0/query?outFields=*&where=1%3D1&f=geojson')

Reading layer `OGRGeoJSON' from data source 
  `https://atoll.floridamarine.org/arcgis/rest/services/FWC_GIS/OpenData_Shoreline/MapServer/0/query?outFields=*&where=1%3D1&f=geojson' 
  using driver `GeoJSON'
Simple feature collection with 15470 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -87.63477 ymin: 24.52114 xmax: -80.03118 ymax: 31.00091
Geodetic CRS:  WGS 84

flshore <- st_transform(flshore, crs = 6443)

dattbdbasins <- dattbdbasins %>% 
  st_intersection(flshore)

dattbindex <- dattbindex %>% 
  st_intersection(flshore)

Calculate (1) the total land area of each drainage basin in the watershed and (2) the total area of underserved tracts within each drainage basin. Units are in feet for this projected coordinate system.

# get only underserved tracts
dattbunder <- dattbindex %>%
  filter(underserved == "Yes")

# total basin area
basinareas <- dattbdbasins %>%
  mutate(area_ft = st_area(dattbdbasins)) %>%
  group_by(WBID) %>%
  summarise(db_area_ft = sum(area_ft))

# underserved area within basins
basin_int <- st_intersection(basinareas, dattbunder)

basin_underareas <- basin_int %>%
  mutate(area_int_ft = st_area(basin_int)) %>%
  group_by(WBID) %>%
  summarise(under_area_ft = sum(area_int_ft)) %>%
  as.data.frame()

# join area estimates
areasjoined <- left_join(basinareas, basin_underareas, by = 'WBID') %>%
  rowwise() %>%
  select(matches('^WBID|^db|^under')) %>%
  as.data.frame()

dattbdbasins_under <- left_join(dattbdbasins, areasjoined, by = 'WBID')

Calculate the proportion of each drainage basin containing underserved communities.

dattbdbasins_under <- dattbdbasins_under %>%
  mutate(under_area_ft = ifelse(is.na(under_area_ft), 0, under_area_ft)) %>%
  mutate(pct_under = under_area_ft/db_area_ft * 100) %>%
  mutate(pct_under = as.numeric(pct_under))

Create a map showing priority drainage basins based on the presence of underserved communities. Hover over the drainage basins to view the name of the water body assessment unit.

mapviewOptions("basemaps.color.shuffle" = FALSE)

mapview(dattbdbasins_under, zcol = "pct_under", col.regions = brewer.pal(6, "Reds"), at = seq(0, 100, 20), layer.name = "Underserved Tracts (% of DB)", label = "WATERBODY_NAME")